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ABSTRACT 

We develop a simple dynamical model for the evolution of gas in the centres of barred spiral 
galaxies, using the Milky Way’s Central Molecular Zone (CMZ, i.e., the central few hundred 
pc) as a case study. We show that, in the presence of a galactic bar, gas in a disc in the 
central regions of a galaxy will be driven inwards by angular momentum transport induced by 
acoustic instabilities within the bar’s inner Lindblad resonance. This transport process drives 
turbulence within the gas that temporarily keeps it strongly gravitationally stable and prevents 
the onset of rapid star formation. However, at some point the rotation curve must transition 
from approximately flat to approximately solid body, and the resulting reduction in shear 
reduces the transport rates and causes gas to build up, eventually producing a gravitationally- 
unstable region that is subject to rapid and violent star formation. For the observed rotation 
curve of the Milky Way, the accumulation happens ~ 100 pc from the centre of the Galaxy, in 
good agreement with the observed location of gas clouds and young star clusters in the CMZ. 
The characteristic timescale for gas accumulation and star formation is of order 10 — 20 Myr. 
We argue that similar phenomena should be ubiquitous in other barred spiral galaxies. 

Key words: galaxies: evolution — Galaxy: centre — Galaxy: evolution — ISM: kinematics 
and dynamics — stars: formation 


1 INTRODUCTION 

Star formation is one of the most important unsolved problems in 
contemporary astrophysics, and one of the central uncertainties is 
the extent to which the rate of star formation per unit mass of inter¬ 
stellar gas is influenced by galactic-scale processes. (See Krumholz 
2014, Dobbs et al. 2014, and Padoan et al. 2014 for recent re¬ 
views.) On one hand, a number of observations appear to favor a 
“bottom-up” view in which star formation is a purely local pro¬ 
cess that does not depend on galactic characteristics. These include 
the apparent insensitivity of the star formation rates in galaxies to 
galactic parameters such as Toomre (1964) Q (e.g., Leroy et al. 
2008, 2013) or to the presence and structure of spiral arms (Willett 
et al. 2015), the minimal variation in molecular gas depletion time 
(fdep = Mgas/SFR) measured on ~ kpc scales in nearby galaxies 
(e.g., Bigiel et al. 2008; Schruba et al. 2011; Leroy et al. 2013), 
and the fact that all star-forming systems, on scales from individual 
clouds to entire starburst galaxies, appear to turn their mass into 
stars at a nearly constant rate per gas free-fall time (Krumholz & 
Tan 2007; Krumholz, Dekel & McKee 2012), possibly with a small 
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secondary variation based on the Mach number of the turbulence 
(Federrath 2013). In this view, the global rate of star formation in 
galaxies is mostly a matter of adding up individual star-forming 
clouds, whose behavior is independent of their galactic environ¬ 
ment. Models based on this premise have been advanced by a num¬ 
ber of authors (e.g., Krumholz & McKee 2005; Krumholz, McKee 
& Tumlinson 2009; Lada et al. 2012; Renaud, Kraljic & Bournaud 
2012; Krumholz 2013; Federrath 2013; Salim, Federrath & Kewley 
2015). 

On the other hand, there is a broad class of theoretical models 
predicts that galaxies’ star formation rates should explicitly depend 
on galactic structure, at least in some galaxies (e.g., Thompson, 
Quataert & Murray 2005; Ostriker, McKee & Leroy 2010; Ostriker 
& Shetty 2011; Shetty & Ostriker 2012; Hopkins, Quataert & Mur¬ 
ray 2011; Faucher-Giguere, Quataert & Hopkins 2013), and a num¬ 
ber of recent observations using higher resolution data have pointed 
to a somewhat greater role for galactic dynamics. These include 
observations showing that the properties of molecular clouds vary 
systematically with environment both between and within galax¬ 
ies (e.g., Hughes et al. 2013; Colombo et al. 2014), and evidence 
that the star formation rate per unit molecular mass measured on 
small scales is not independent of the large-scale rate of shear in a 
galactic disc (Meidt et al. 2013; Suwannajak, Tan & Leroy 2014). 


© 2015 RAS 


2 Krumholz & Kruijssen 


A particularly promising avenue for exploring the role of 
galaxy-scale dynamics in regulating star formation is to study the 
centres of galaxies. These are the regions where shear and simi¬ 
lar effects arising from the shape of the galactic potential should 
have their largest effects (e.g., Kruijssen et al. 2014; Kruijssen, 
Dale & Longmore 2015). It is also the location where other envi¬ 
ronmental effects, such as high external pressures (e.g., Rathbome 
et al. 2014b) and high X-ray and cosmic ray fluxes (e.g., Meijerink 
& Spaans 2005; Papadopoulos 2010; Meijerink et al. 2011; Clark 
et al. 2013; Kruijssen et al. 2014), should be at their strongest. 
Indeed, galactic centres show a number of interesting deviations 
from the star formation behavior seen at larger galactic radii. Some 
galaxies show enhanced star formation per unit molecular mass in 
their nuclei, others depressed rates of star formation (Saintonge 
et al. 2012; Leroy et al. 2013; Longmore et al. 2013a). Galactic 
centre molecular clouds show systematic differences in their prop¬ 
erties from disc clouds (e.g., Kruijssen & Longmore 2013; Leroy 
et al. 2014; Rathborne et al. 2014b, 2015; Bally et al. 2014). 

The best-studied of galactic centres is the Central Molecular 
Zone (CMZ) of the Milky Way, which exhibits a number of inter¬ 
esting features. Much of the dense gas in this region appears to be 
collected into a stream, or a partially filled ring, of molecular mate¬ 
rial ~ 100 pc from the Galactic centre (Molinari et al. 2011; Kruijs¬ 
sen, Dale & Longmore 2015). Within this ring are a series of clouds 
whose star forming activity ranges from far smaller than one would 
expect based on their mass and density (e.g., GO.253+0.016, also 
known as “The Brick”; Longmore et al. 2012; Kauffmann, Pillai & 
Zhang 2013; Rathborne et al. 2014a; Mills et al. 2015) to some of 
the most actively star-forming sites in the Local Group (Sgr A, Sgr 
B2, and Sgr C; Yusef-Zadeh et al. 2008, 2009). These clouds may 
well represent an evolutionary sequence (Longmore et al. 2013b; 
Kruijssen, Dale & Longmore 2015). There is also significant evi¬ 
dence that star formation in the CMZ is episodic. Hints at previ¬ 
ous starbursts include the presence of large off-plane bubbles vis¬ 
ible in radio (Sofue & Handa 1984), infrared (Bland-Hawthom & 
Cohen 2003), and gamma-rays (Su, Slatyer & Finkbeiner 2010), 
along with direct counts of young stellar objects (Yusef-Zadeh et al. 
2009). 

In this paper we propose a simple model for the global dynam¬ 
ics and star formation behaviour of gas in the Milky Way’s CMZ, 
and by extension the central regions of other barred spiral galaxies. 
Qualitatively, our model is that gas is channelled through the disc 
of the Milky Way, along the Galactic bar, and into the CMZ. Once 
there, it settles into a disc that is subject to acoustic instabilities 
driven by the bar. These instabilities drive both turbulence and an¬ 
gular momentum transport in the gas, simultaneously causing it to 
flow inward and increase in velocity dispersion. As a result, the gas 
is extremely turbulent, and its large velocity dispersion renders it 
highly stable against gravitational collapse or star formation (typi¬ 
cally Q ~ 100, as we show below). However, near the radius where 
the rotation curve of the Galaxy turns over from flat to solid body, 
the resulting reduction in shear suppresses transport and turbulent 
driving, leading gas to accumulate rather than moving further in¬ 
ward. This accumulation eventually builds up a ring of material 
which goes gravitationally unstable and begins vigorous star for¬ 
mation. If feedback from this burst of star formation is sufficiently 
strong, the ring will be disrupted, and the cycle will begin again. 

Our plan for the remainder of this paper is as follows. In Sec¬ 
tion 2 we introduce our physical model for the gas in the CMZ, 
and in Section 3 we describe how we simulate the evolution of this 
model. In Section 4 we present the results of our simulations, and 


we discuss the implications of those results in Section 5. We sum¬ 
marize and conclude in Section 6. 


2 MODEL 

Our goal in this section is to develop a quantitative model corre¬ 
sponding to the qualitative scenario described in the previous sec¬ 
tion. We will model the gas near a galactic centre as an axisym- 
metric thin disc subject to non-axisymmetric perturbations, within 
which mass and angular momentum are transported via instabilities 
and energy is lost due to the dissipation of supersonic turbulence. 
We will simulate this system using the VADER code of Krumholz & 
Forbes (2015). In the remainder of this section we detail the physi¬ 
cal ingredients to our model, and in the following one we describe 
the numerical setup we use to simulate it. Our simulations will fo¬ 
cus on the case of the Milky Way CMZ, since it is the system for 
which we have the best measurements of the small-scale rotation 
curve. 


2.1 The Galactic Potential and the Inner Lindblad 
Resonance 


The first step in this modeling is to derive the potential in which the 
gas orbits. To do so, we make use of two data sets for the Milky 
Way’s CMZ. Our primary data set, which we use for all numerical 
calculations, is the enclosed mass M r versus galactocentric radius r 
measured by Launhardt, Zylka & Mezger (2002) at r = 0.63 — 488 
pc, which we use to derive a rotation speed = yjGM r /r. Our 
secondary data set, which we will not use for numerical compu¬ 
tations but that we include here because the Launhardt, Zylka & 
Mezger (2002) data do not go out far enough to reach the inner 
Lindblad resonance (ILR; e.g., Binney & Tremaine 1987), is the 
rotation curve compiled by Bhattacharjee, Chaudhury & Kundu 
(2014) from r = 190 — 1.9 x 10 5 pc. 1 We can interpolate these 
tabulated data to generate a rotation curve v$, versus r, but this re¬ 
quires some care. The dynamical evolution of gas in this potential 
depends not just on the rotation curve but on its gradient, 


P 


d In Vtf, 

dlnr 


( 1 ) 


This quantity is important because the dimensionless rate of shear 
is 1 — /?. For this reason we interpolate using basis splines (B- 
splines), following the method of Gans & Gill (1984), as imple¬ 
mented in the VADER code by Krumholz & Forbes (2015). The 
B-spline fit ensures continuity of some number of derivatives. For 
the Launhardt, Zylka & Mezger (2002) data set we use 6th order 
B-splines with 15 breakpoints, while for the Bhattacharjee, Chaud¬ 
hury & Kundu (2014) data set, since it contains fewer data points 
within 6 kpc (where we truncate the fit), we use 3rd order B-splines 
with 6 breakpoints. The top panel of Figure 1 shows the data and 
our fits to them. We see that there is some tension between the two 
data sets, but that there is rough qualitative agreement in the range 
of radii where they overlap. The Figure also illustrates why the B- 
spline fit is critical: a simple linear fit produces rates of shear with 
unphysical sharp features. 


1 Bhattacharjee. Chaudhury & Kundu (2014) provide three possible fits, 
corresponding to three choices for the imperfectly known Galactocentric 
radius of the Solar Circle. We use their data set corresponding to 8.0 kpc for 
this value, but the results are qualitatively identical for other choices. 
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Figure 1. Rotation curve for the inner Galaxy. Top panel: rotation speed 
versus Galactocentric radius r using the data set of Launhardt, Zylka 
& Mezger (2002, blue circles and thick solid line) and of Bhattacharjee, 
Chaudhury & Kundu (2014, green squares and thin solid line). For the 
versus r plots, the points indicate the observational estimate, while the lines 
are B-spline fits to the data (see text for details). The dashed blue line shows 
the dimensionless shear 1 — ft = 1 — din /din r derived from the B 
spline fit to the Launhardt, Zylka & Mezger data; for comparison, the dotted 
blue line shows 1 — ft computed using a simple linear fit. Bottom panel: 
angular frequencies versus radius. Solid lines show Qq — re/2. where CIq 
is the orbital frequency and re is the epicyclic frequency, for the fits to the 
Launhardt, Zylka & Mezger ( thick blue line) and Bhattacharjee, Chaudhury 
& Kundu ( thin greeti line) data sets. The thin black dashed line shows the 
pattern speed of the Galactic bar; the intersection of this line with Oo — re/2 
is location of the ILR. 


The lower panel of Figure 1 shows, for our fits to both data 
sets, the frequency f 1 — k/2, where fl = v^/r is the angular 
velocity of the orbit and re = \/2(l + /3)0 is the epicyclic fre¬ 
quency. For comparison, the black dashed line shows the pattern 
speed fibar ~ 0.06 Myr” 1 of the Galactic bar (Debattista, Gerhard 
& Sevenster 2002; Wang et al. 2012; Antoja et al. 2014); the point 
where Q, — k/ 2 = S2bar marks the location of the ILR. We see that 
the ILR is located at ~ 1 kpc. The entire region we model will be 
within this radius. 

Finally, note that, on top of this cylindrically-symmetric po¬ 
tential, there is an additional non-axisymmetric component to the 
potential due the presence of a bar. We incorporate this effect 
into our model below by considering how the presence of a non- 
axisymmetric perturbation affects this disc, and so we defer any 
further considering of the bar’s effects for the moment. 


2.2 Gas Evolution 


We approximate that gas in the CMZ orbits in an azimuthally- 
symmetric thin disc, which is characterized by a surface density 
E, a non-thermal velocity dispersion a nt , and a thermal veloc¬ 
ity dispersion crth- Both E and <j n t are functions of r, but we set 
a t h = 0.5 km s -1 independent of position. This is the sound speed 
of fully molecular gas at a temperature of 70 K, in the middle of the 
kinetic temperature range for CMZ gas that Ao et al. (2013, also 
see Ginsburg et al. 2015) derive by analyzing H 2 CO line emission. 
The exact choice will have little effect on our results in any event, 
because observations show, and the model we describe below also 
predicts, that a nt cr t h almost everywhere in the disc. 

The gas evolves following the standard equations of mass and 
energy conservation for such a disc (e.g., Krumholz & Burkert 
2010 ), 


-E+ - — 
dt r dr 

d — 1 fl , ,_ Id 

dt E+ r!fr\ rVr ^ + ~ r <9r 


(rv r E) 
VT 3 

2nr 2 ) 


0 (2) 
-E ra d 5 (3) 


where E = E \ip + r \/2 + (3/2) + of h )] is the total en¬ 

ergy of the gas per unit area, tp is the gravitational potential, 
.E/ad is the rate of energy gain or loss due to radiative effects, 
P = E (ffj t + <7t h ) is the total vertically-integrated pressure, v r 
is the radial velocity, and T is the turbulent torque. The latter two 
quantities are related via angular momentum conservation, 


dT/dr 

27rrEu0(l + ft) 


(4) 


We parameterize the viscosity though the Shakura & Sunyaev 
(1973) a-disc model, including the generalization proposed by Shu 
(1992): 


T = —2ttr 2 aP (1 — /3), 


(5) 


where a is a dimensionless number characterizing the strength of 
angular momentum transport. We defer a discussion of the choice 
of a to the next section. 

It is worth pausing to remark on the physical processes that 
are responsible for controlling the turbulence in the gas that are 
captured in equation 3, since these will be crucial to understanding 
our results. The non-thermal velocity dispersion a nt is a component 
of the total energy per unit mass E/E. Since all the other compo¬ 
nents depend on quantities that are fixed at each radius (ip, v</,, and 
crth), any change in the energy E at a given location is expressed 
as a change cr nt . These changes can be driven by several processes, 
each corresponding to a different term in equation 3. The easiest 
to understand is the decay of supersonic turbulence via radiative 
shocks, which is captured by the term E ra d- This process causes 
(T nt to decrease at every point in the disc. 

The countervailing processes are represented by the 
terms on the left-hand side of equation 3. The term 
— (l/r)(d/dr)(rv c f, / T/2nr 2 ) describes the work done by 
torques. Physically, in our viscous disc model angular momentum 
transport occurs because each ring of material has a higher angular 
velocity than the ring immediately outside it, and as the fluid 
elements move past one another they exert torques. This torque 
does work, the effect of which is to transport energy outward. 
Finally, the term (1 / r)(d / dr)[rv r {E + P )] combines advection 
of energy by bulk motion of gas and work done by pressure 
forces. Since E )$> P for a thin disc, advection is by far the more 
important of these two processes. Advection tends to increase <r n t 
because, although fluid elements retain constant specific energy 
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E /E as they advect, the gravitational potential ip and rotation 
curve -U 0 are such that ip + v^/2 decrease inward. Thus the specific 
orbital plus gravitational energy of a fluid element goes down as 
it is transported inward, and by energy conservation this leads the 
specific turbulent energy (3/2)<r 2 t to increase. Thus the velocity 
dispersion in our disc is dictated by a competition between the 
decay of turbulence in radiative shocks, which decreases <j n t, and 
the combined effects of turbulent torques and advection, which 
tend to increase it by trading off gravitational potential energy 
against kinetic energy. 


2.3 Energy Dissipation Processes 

The next step in our model is to evaluate E la a. In the absence of 
any forcing, the non-thermal velocity dispersion of the gas will de¬ 
cay with time. Numerical simulations indicate that the timescale 
for turbulence to decay is of order the crossing time of the flow 
(Stone, Ostriker & Gammie 1998; Mac Low 1999; Ostriker, Stone 
& Gammie 2001; Lemaster & Stone 2009). In the context of a 
disc, this timescale should be computed relative to the gas scale 
height, which we obtain using the model of Ostriker, McKee & 
Leroy (2010) for the vertical gas structure. In this model, the scale 
height H g of the gas is implicitly given by 

2n( d Gp t ZH 2 + ^GE 2 H g = P, (6) 

where p, is the stellar (plus dark matter) mass density within a 
distance H g of the gas midplane (assumed to vary negligibly, which 
is true as long as the stellar scale height greatly exceeds H g ) and 
C,d 0.33 is a dimensionless constant whose exact value depends 
on the relative importance of gas self-gravity versus stellar gravity, 
but whose value remains within 5% of 0.33 in all cases. Intuitively, 
equation 6 simply asserts that, in hydrostatic balance, the vertically- 
integrated pressure must balance the gravitational forces attempting 
to compress the gas, which come from the gravitational potential 
produced by stars and dark matter (the first term on the left hand 
side) and the self-gravity of the gas (the second term). 

To obtain the stellar density, we note that, if the stellar mass 
responsible for producing the rotation velocity were arranged 
spherically, we would have 

P* = /shape (1 + 2/3) ^Q r 2 (7) 

with /shape = 1. A flattened stellar distribution has /shape > 1. 
We do not precisely know the full three-dimensional distribution of 
mass at the Galactic Centre, but /shape can be constrained some¬ 
what by observations. Rodriguez-Fernandez & Combes (2008), 
based on 2MASS star counts, and Molinari et al. (2011), based on 
the shape of the 100 pc gas structure in the CMZ, both estimate 
/shape ~ 2. More recently, Kruijssen, Dale & Longmore (2015) 
fit the potential based on the kinematics of gas in the CMZ 2 , and 
obtained / s ha P e ~ 2.5. We will adopt this values as our fiducial 
choice, but also study how varying it affects our results. 

Given a choice for /shape, at any point in the disc with known 
E, P, and v^, equation 6 is a quadratic with known coefficients that 
we can solve to obtain H g . We then set the energy dissipation rate 
to 


2 Kruijssen, Dale & Longmore (2015) express their result in terms of the 
vertical shape of the potential, expressed via the axial ratio of the equipo- 
tential surfaces q,\, (Binney & Tremaine 1987, section 2.2.2). Their favored 
value is equivalent to /shape = 2.5. 


E Tad — Tj 


Ea 2 


TTg/Cnt ’ 


(8) 


where 77 is a constant of order unity, for which we take a fiducial 
value of 1.5, calibrated from simulations (Stone, Ostriker & Gam¬ 
mie 1998). This amounts to setting the rate of energy loss equal to 
the current kinetic energy divided by one gas scale-height crossing 
time. 


2.4 Angular Momentum Transport Processes 


The final important element in our model is angular momentum 
transport, as parameterized by the dimensionless value a. To es¬ 
timate this quantity, consider a disc of molecular gas in a barred 
galaxy, such as that in the CMZ. The stellar bar will perturb the 
gas disc with an m = 2 mode, and we wish to consider the sta¬ 
bility of the disc against these perturbations. Stability analyses 
such as this have been performed by a large number of authors 
(e.g., Goldreich & Lynden-Bell 1965; Lau & Bertin 1978; Toomre 
1981; Bertin et al. 1989; Montenegro, Yuan & Elmegreen 1999). 
Bertin et al. (1989) obtained the general dispersion relation for non- 
axisymmetric self-gravitating thin discs, and Montenegro, Yuan & 
Elmegreen (1999) pointed out that it admits two classes of instabil¬ 
ity: gravitational instabilities, which occur for axisymmetric per¬ 
turbations and for non-axisymmetric ones close to co-rotation, and 
acoustic instabilities, which arise strictly for non-axisymmetric per¬ 
turbations inside the ILR or outside the outer Lindblad resonance 
(OLR). Acoustic instabilities are so named because the destabiliza¬ 
tion is driven by pressure rather than gravity; pressure causes the 
apocentres of perturbed gas orbits inside the ILR to align, lead¬ 
ing to a growing mode. Indeed, gravity acts as a stabilizer rather 
than a destabilizer for acoustic instabilities, since it “de-tunes” the 
pressure-driven alignment. 

Formally, the dispersion relation for waves in a self- 
gravitating thin disc is (Montenegro, Yuan & Elmegreen 1999, their 
Equation 10) 


Cf 

4 


V ■ 


1 - v 2 


fj~ 2 + J 2 /(l — v 2 ) 


(9) 


where v = (u> — mQ.)/K is the dimensionless frequency of the 
perturbation, uj is the frequency, m is the azimuthal wavenumber 
of the perturbation, and i) is a dimensionless inverse wavenumber, 
which must be positive and real. The instability is controlled by the 
two parameters Q and J, where 


Q = - 

V ttGE 

is the usual Toomre (1964) parameter, and 


J = 


Vn 

&crit 


( 10 ) 


( 11 ) 


with Tj = (1 — /3)(2mfi/ nr) 2 and fe crit = k 2 /27tGE, is a di¬ 
mensionless measure of the strength of shear in the disc. Note that 
J = 0 for /3 = 1, corresponding to solid body rotation and thus 
zero shear. The disc is unstable to the growth of perturbations with 
azimuthal wavenumber m at a given combination of Q and J if 
there exists a positive real value of f) such v has a non-zero imagi¬ 
nary part. 


Equation (9) can be solved for v as 

1 ( Q 2 1 

V 




( 12 ) 


with 
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D = 




l 

v 


A 7-2 -2 

— 4J rj 


(13) 


However, be warned that this re-arrangement introduces some 
spurious solutions, as is often the case with rational equations. 
Nonetheless, this re-arrangement is useful, because it clearly shows 
how to analyze a disc to determine whether it is either gravitation¬ 
ally or acoustically unstable. 

The two regimes of instability correspond to different signs of 
the discriminant D. The gravitational instability case arises when 
D > 0, so v 2 is a real number, but v 2 < 0 so that v is a pure 
imaginary. In this case the transition from stability to instability 
occurs in the vicinity of v = 0 , which corresponds to unstable 
waves being close to co-rotation with the disc (w ~ mfi). For 
m = J = 0, corresponding to axisymmetric perturbations, the 
dispersion relation (equation 9) reduces to 


= 1 + 

4?7 


(14) 


which has a minimum of v 2 = 1 — Q 2 at rj = Q 2 /2. This 
minimum is negative, indicating instability, if Q < 1, the usual 
Toomre (1964) stability condition. For non-axisymmetric perturba¬ 
tions, to 7 ^ 0 , but in real astrophysical systems we generally have 
J < 1 (see Bertin et al. 1989 for more discussion of this point), 
so this condition is modified only slightly. Acoustic instability cor¬ 
responds to the regime where D is itself negative, and thus v 2 and 
v both have non-zero imaginary parts. In this case the real part of 
v is either < — 1 or > 1 , and thus in the portion of the disc that 
lies either inside the ILR or outside the OLR with respect to the 
perturbation. 

For practical numerical purposes, we can evaluate the suscep¬ 
tibility of a disc with specified values of E, a, and rotation curve 
D 0 to acoustic instability via the following procedure. First, we find 
the value of fj that minimizes D\ the global minimum of D, if it 
is less than zero, corresponds to the fastest growing acoustic mode. 
Local minima of D occur at values of fj that satisfy dD/dr) = 0, 
and evaluating the derivative gives 


16 J 2 if — 8fj 2 + 6Q 2 fj — Q 4 = 0. 


(15) 


We obtain the roots of this polynomial via standard numerical tech¬ 
niques (Galassi et al. 2009), discard roots where R e(fj) < 0 or 
Im( 7 )) 0 (since only positive real roots are physically allowed), 

and check the remaining roots to see which one minimizes D. If D 
is negative at its minimum, the region is acoustically unstable, and 
we can obtain the growth rate of the most unstable mode simply 
by plugging this minimum into equation 9 and finding which of the 
four possible values of v has the largest imaginary part. 

Formally, we define the growth time of the fastest growing 
mode as 

t growth — TV j x, ; ( 16 ) 

/tmax[Im(i/)| 

where in evaluating possible values of v we consider both gravita¬ 
tional and acoustic instabilities. To investigate the stability of the 
Milky Way’s CMZ, we use the Launhardt, Zylka & Mezger (2002) 
potential, again treating the potential as cylindrically-symmetric, 
and regarding the much smaller non-axisymmetric component pro¬ 
vided by the bar as a perturbation. For this potential, we compute 
the stability of the disc against both m — 0 gravitational instabil¬ 
ities and to = 2 acoustic perturbations driven by the Galactic bar, 
using a range of sample values of E and a. We show the results 
of this computation in Figure 2. We find that, for the parameter 
choices shown, which span the plausible range based on observa- 



r [pc] 

Figure 2. Growth timescale f grow th for the fastest-growing acoustic in¬ 
stability mode at the Galactic Centre, normalized to the local orbital pe¬ 
riod f or b = 2n/fl, as a function of Galactocentric radius r. The growth 
timescales shown have been computed for values of the surface density 
E = 10 and 100 Mq pc -2 and velocity dispersion <r = 20 and 80 km 
s -1 , independent of position, and for m = 2 modes. The rotation curve 
used in the computation is produced by B-spline interpolation of the Laun¬ 
hardt, Zylka & Mezger (2002) measurements following the procedure out¬ 
lined in Section 2.1, as shown by the thick blue solid line in the upper panel 
of Figure 1. 


tions, the disc is unstable to acoustic modes at all radii; it is not 
unstable to gravitational modes anywhere, although, as we will see 
below, it will naturally evolve into a gravitationally unstable state. 
The growth times of the acoustic modes are generally compara¬ 
ble to the orbital period, indicating that the instability should grow 
efficiently. The exception is just inside ~ 100 pc, where the rota¬ 
tion curve approaches solid body and the shear diminishes (c.f. Fig¬ 
ure 1). This will prove to be important below. 

The presence of an instability will certainly produce turbu¬ 
lence and angular momentum transport. Unfortunately, we cannot 
easily compute the exact transport rate once the instability reaches 
full non-linear saturation. However, it seems likely that they will 
be high. Simulations of both pure gas discs and gas plus stellar 
discs show that gravitational instabilities tend to produce transport 
rates corresponding to a dimensionless viscosity a ~ 1 (e.g., Bour- 
naud, Elmegreen & Elmegreen 2007; Kratter et al. 2010; Hopkins 
& Quataert 2010; Ceverino et al. 2015). Acoustic instabilities, un¬ 
like gravitational ones, cannot self-stabilize by producing turbu¬ 
lence and thus driving up the value of Q (e.g., Krumholz & Burkert 
2010), so the instability seems likely to be at least as strong. It is 
likely that the growth of the instability is ultimately limited by the 
dissipation of energy in spiral shocks. 

Given the lack of a well-motivated estimate for the transport 
rate due to acoustic instability, we choose to parameterize it by 

a = min(aoe 1 ~ tErowth/,t ° rb , 1) (17) 

where t or b = 2n/Cl is the local orbital period. The parameter 
ao simply normalizes the rate of angular momentum transport at 
points where the instability grows on a timescale of a single orbital 
period, while taking the minimum ensures that, even in the most 
unstable regions, we do not exceed a = 1. We will adopt a fidu¬ 
cial value ao = 1 , but consider below how varying this parameter 
might change the results. 


© 2015 RAS, MNRAS 000, 1-19 










6 Krumholz & Kruijssen 

2.5 Summary of the Model 

We now have in place a fully-specified model of the evolution of 
gas in the Milky Way’s CMZ. The overall evolution is described 
by equation 2 and equation 3, which specify mass and energy con¬ 
servation. These equations depend on the rotation curve v# and po¬ 
tential <j>, the rate of energy dissipation _E ra d, and the dimension¬ 
less rate of angular momentum transport a. For the rotation curve 
and potential we use the B-spline fit to the Launhardt, Zylka & 
Mezger (2002) rotation curve shown in Figure 1. For the rate of 
radiative energy loss we use the value given by equation 8, which 
comes from our model of turbulent dissipation. Finally, for angular 
momentum transport we use our calculation of acoustic and grav¬ 
itational instabilities, as parameterized by equation 17. This set of 
equations, once we choose values for the various parameters ap¬ 
pearing in them, fully specifies the problem. 


3 SIMULATION METHOD AND SETUP 

We simulate the disc model described in the previous section us¬ 
ing the VADER code described by Krumholz & Forbes (2015). 3 
VADER solves the equations using a fully-implicit method that con¬ 
serves mass and energy to machine precision. We use VADER’s 
backwards-Euler solver, because the Crank-Nicolson one produces 
undesirable numerical oscillations. All our simulations use a com¬ 
putational grid of 512 cells, spaced logarithmically from an inner 
edge at r = 10 pc to an outer edge at r = 450 pc. Although 
the region inside the ILR extends to ~ 1 kpc, and the bar goes 
out even further, we choose 450 pc for the outer edge because 
it is close to the outermost point included in Launhardt, Zylka & 
Mezger (2002)’s estimate of the Galactic rotation curve, and thus 
avoids the need to interpolate between this rotation curve and that 
of Bhattacharjee, Chaudhury & Kundu (2014). We have verified 
that choosing different values that are still of order hundreds of pc 
does not qualitatively alter the results. 

The equations describing the evolution of the disc require two 
boundary conditions at each end, describing the rate at which mass 
and enthalpy are advected into or out of the computational domain. 
At the inner boundary, we prescribe our boundary by requiring that 
the viscous torque go to zero. This implicitly sets the mass flux 
(see Krumholz & Forbes 2015 for details), and forces it to be out of 
the computational domain. We set the specific enthalpy at the inner 
boundary equal to the initial specific enthalpy in the simulation (see 
below), but this has no practical effect since no material ever enters 
the computational domain at the inner edge. 

The choice of outer boundary condition is less trivial. We 
choose to specify the boundary condition as a fixed inward mass 
flux, but we do not have good observational estimates of the rate at 
which material is transported through the outer parts of the Milky 
Way’s disc to enter the CMZ, or of the velocity dispersion of this 
material as it arrives. Theoretical estimates suggest that the inward 
transport rate through the disc should be of order Mj n = 1 Mg 
yr _1 (Krumholz & Burkett 2010; Forbes, Krumholz & Burkert 

3 We make one minor modification to the system as described: we adopt 
a minimum viscosity a = 10 -3 . We do this both because very small val¬ 
ues of a cause difficulties in convergence that slow the calculations, and 
because we expect a to be at least this large as a result of thermal and 
magneto-rotational instability (e.g. Piontek & Ostriker 2004, 2007), even 
in regions that are acoustically and gravitationally stable. The value of the 
floor parameter does not affect the results. 


2012; Forbes et al. 2014; Cacciato, Dekel & Genel 2012), and that, 
once gas reaches the bar, it should be further transported inward 
along the bar to be deposited in the CMZ (Binney et al. 1991; Kor- 
mendy & Kennicutt 2004), where it settles into the disk-like struc¬ 
ture that we seek to model (Sormani, Binney & Magorrian 2015). 
We therefore adopt M ln = 1 Mg yr _1 as a fiducial value, but be¬ 
low we will explore how robust our results are against variations in 
this choice. 

Observations indicate that the velocity dispersion of the dens¬ 
est molecular gas near the galactic centre is of order tens of km 
s _1 (e.g., Walsh et al. 2011; Purcell et al. 2012, and Table 1 of 
Kruijssen et al. 2014 4 ). The velocity dispersion of somewhat lower 
density gas is much less constrained, but observations elsewhere 
in the Galaxy (Walsh, Myers & Burton 2004; Andre et al. 2007; 
Kirk, Johnstone & Tafalla 2007; Rosolowsky et al. 2008), as well 
as theoretical expectations (e.g. Padoan et al. 2001; Offner et al. 
2008; Offner, Hansen & Krumholz 2009), suggest that it should 
be larger. We adopt crj n = 40 km s _1 as our fiducial velocity dis¬ 
persion, but, as with the mass inflow rate, we explore below how 
varying ai n changes our results. 

As initial conditions for our simulations, we choose to place 
a small amount of material into the computational domain. Specif¬ 
ically, we set the initial surface density to S = 1 Mg pc -2 and 
the initial velocity dispersion to cr = 40 km s _1 everywhere. The 
results are quite insensitive to these values as long as the mass and 
energy of the initial gas is small enough that the external inflow 
rapidly swamps the mass and internal energy of the material that 
was present at the beginning. For our fiducial parameter choices, 
externally-input material dominates after 0.6 Myr. 

All the source code for the simulations described in 
this section, and for the analysis presented in the subse¬ 
quent sections, is publicly available. The run and analysis 
scripts can be obtained from https://bitbucket.org/ 
krumholz/cmzdisk/ (hash b98da00), and VADER code is 
available at https://bitbucket.org/krumholz/vader 
(hash 5a96350). 

4 RESULTS 
4.1 Fiducial Case 

We first use the numerical method described in the previous section 
to simulate a fiducial case, for which the free parameters are as 
shown in Table 1. Figure 3 shows the results of this simulation over 
50 Myr. Qualitatively, we see that mass enters the computational 
domain from the outer edge and propagates inward through the disc 
as time passes. At early times the inflow rate is highest at the outer 
edge of the disc (where it is forced by our boundary conditions to 
be 1 Mg yr _1 ), and then decreases inward. The velocity dispersion 
rises as the gravitational potential energy of this incoming material 
is converted into turbulent motion, at a rate that is initially too high 
for the decay of turbulence to counter it. This increase in velocity 
dispersion keeps Q high at radii above ~ 100 pc. By ~ 30 Myr of 
evolution the disc outside ~ 100 pc has settled into a steady state 
whereby the inflow rate produced by acoustic instability is the same 
at all radii, and the velocity dispersion is kept constant by a balance 
between turbulent dissipation and the conversion of gravitational to 

4 These authors generally report either the FWHM or the ID velocity dis¬ 
persion, while our a is the 3D velocity dispersion, which is larger by a factor 
of v/3. 
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Table 1. Parameter values for the fiducial case. 


Parameter 

Value 

Meaning 

Q!0 

1.0 

Dimensionless viscosity 

Min 

1.0 M 0 yr" 1 

External accretion rate 


40 km s -1 

Incoming 3D velocity dispersion 

/shape 

2.5 

Potential shape parameter 

^th 

0.5 km s -1 

Thermal velocity dispersion 


1.5 

Turbulence decay rate 


turbulent energy as mass accretes down the potential well. In this 
equilibrium, Q remains of order 100, so that the gas is extremely 
stable against collapse or star formation. 

The results change qualitatively inside ~ 100 pc, as the wave 
of incoming material stops propagating inward and instead piles up 
in a ring. The inward mass flow rate also drops sharply as acous¬ 
tic instability shuts off. At later times gravitational instability takes 
over, but this does not prevent a large pileup of gas. Within the ac¬ 
cumulating ring of material, the velocity dispersion begins to drop, 
reaching the thermal floor value of 0.5 km s -1 . As a result the Q 
value in the ring begins to drop, first falling below unity at some 
point within the disc roughly 10 — 15 Myr after the simulation 
starts. In Figure 4 we show the total amounts of mass for which 
Q < 1 and Q < 10, and the mass-weighted mean radius of this 
gas, as a function of time. 5 By ~ 50 Myr, several times 10 7 Mg 
of material has accumulated in the unstable region. Thus a mass 
comparable to that of the observed ring-like stream in the CMZ 
(Molinari et al. 201 1) can accumulate in an unstable region after a 
few tens of Myr. 

We can understand the outcome of the simulations by exam¬ 
ining the rotation curve shown in Figure 1. The region where the 
incoming material stalls is precisely where the rotation curve turns 
over from close to flat, as it is in the bulk of the Galaxy, to close to 
solid body, as it is near the Galactic Centre. In this near-solid body 
region, there is little shear, as shown by the dimensionless shear 
1 — /3. This suppresses transport in two ways. First, since acous¬ 
tic instabilities are driven by the presence of shear, the approach to 
solid body rotation weakens the acoustic instability and thus low¬ 
ers a. Second, recalling equation 5, we see that even for fixed a the 
torque is reduced in low-shear regions. This is physically what we 
should expect: for any plausible turbulent transport mechanism, the 
rate of angular momentum transport should be proportional to the 
rate of shear between adjacent rings. If there is no shear, transport 
mechanisms will shut down. 

5 This computation is somewhat subtle, because the regions in which we 
are interested are, at early times, just a few computational cells wide. To 
suppress discreteness noise, we compute the mass in the unstable region via 
the following procedure. We take the simulation output at each time and 
compute Q on our computational grid. We then construct Akima (1970) 
spline approximations to log r versus log £ and log r versus log Q at each 
time step. We use the splines to generate finer grids of 32,768 points, and 
integrate the mass in the Q < 1 and Q < 1(1 regions on that higher resolu¬ 
tion grid. Note that the small oscillations seen in mass versus time shown for 
Q < 1 in Figure 4 are not a result of this procedure, and are robust against 
changes in it. Instead, they are a real feature of the simulation. Oscillations 
occur because as Q approaches unity, gravitational instability turns on as 
a transport mechanism. This leads parts of the disc to undergo a cycle in 
which regions where Q < 1 have active mass transport that raises the ve¬ 
locity dispersion and drives Q above unity so that transport stops. At that 
point, the turbulence begins to decay, so Q eventually drops back below 
unity, turning on transport and starting the cycle again. 



r [pc] 


Figure 3. Results of the simulation with the fiducial parameters. From top 
to bottom, the panels show the surface density £, velocity dispersion <r, 
inward radial velocity —v r , instantaneous inflow rate M; n = — 2irrT,v r , 
Toomre Q parameter, and instability growth time normalized to the orbital 
time (growth /t orb- I n the top four panels, we show results at t = 0 ( thin 
black solid line), t = 10 Myr (thick blue solid line), t = 25 Myr (thick 
green dashed line), and / = 50 Myr (thin red dashed line), as indicated 
in the legend. For the instability growth time in the bottom panel, we show 
the value for the fastest-growing mode, with a solid line indicating that the 
fastest mode is acoustic and a dashed line indicating that it is gravitational. 
Colors and line thicknesses are as in the panels above. 
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Figure 4. Mass in the regions with Q < 1 and Q < 10 (top panel, solid 
blue and dashed green lines, respectively), and mass-weighted mean radius 
of this mass (bottom panel ) versus time in the fiducial simulation. The mean 
radius is plotted only at times when the mass is non-zero. 

A third, less obvious effect of the low-shear region in on the 
energy balance of the disc. Energy transport in the disc occurs 
through two channels or roughly equal importance. Material ad¬ 
verting inward converts potential energy into kinetic energy as it 
does so, raising the local velocity dispersion, and torques within 
the disc transport kinetic energy outward from faster-moving ma¬ 
terial near the centre to slower-moving material further out. Both 
of these processes also tend to shut down in regions of low shear. 
As a result, the rate at which the velocity dispersion of material 
is being driven up by transport processes reaches a minimum in 
the near-solid body region. Since the rate of turbulent dissipation 
does not diminish there, this causes the local velocity dispersion to 
diminish, which in turn reduces the rate of turbulent transport to 
diminish even further. 

As a result of these runaway effects, the low-shear region acts 
like a barrier that stops transport of material inward and energy out¬ 
ward, leading to a buildup of a low velocity dispersion, high surface 
density ring. The accumulation continues until gravitational insta¬ 
bility sets in, at which point it seems likely that star formation will 
ensue and that star formation feedback will increase the velocity 
dispersion and possibly eject material entirely. The time required 
to accumulate enough mass to reach gravitational instability for the 
fiducial parameters is ~ 15 Myr. 

4.2 Dependence on Free Parameters 

We next investigate to what extent the results we have obtained for 
the fiducial case are generic, and to what extent they depend on our 
poorly-known parameters. We summarize the parameters we vary, 
and the values we try, in Table 2. We vary the parameters one at 
a time, leaving all others fixed to the fiducial values indicated in 
Table 1. In all cases we use the same initial conditions and com¬ 
putational grid, and run for 50 Myr. We do not vary the dissipation 
rate coefficient rj, partly because it is reasonably well-calibrated 


Table 2. Variations in parameter values. 


Parameter 

Values 

Meaning 

Q=0 

0.1. 0.5,1.0, 2.0 

Viscosity 

Min 

lO" 0 ' 5 , 1, 10°- 5 Mq yr” 1 

Accretion rate 

Gm 

20, 30, 40, 80 km s -1 

Velocity dispersion 

/shape 

1.0, 2.0, 4.0 

Potential shape 


Fiducial values are indicated in bold. 



Figure 5. Results of the simulations with all parameters set to their fiducial 
values (Table 1 ) except c*o, which varies as indicated in the legend. The top 
row shows the surface density S and the bottom shows Q: the left column 
shows results at t = 15 Myr, and the right at t = 45 Myr. The fiducial case 
is shown as the solid thick black line, in this and all subsequent figures in 
this section. 


from simulations, and partly because it is redundant with /shape- 
We also do not show the results of varying the thermal velocity 
dispersion <r t h, because for plausible values of this parameter its 
effects are negligible. 

We first vary ao, which controls the normalization for the rate 
of angular momentum transport. Figures 5 and 6 show how the re¬ 
sults change as we vary op from its fiducial value of 1.0 in the 
range 0.1 — 2.0. We see that the results for all values of «o -S 0.5 
are qualitatively the same as in the fiducial case. On the other hand, 
for ap = 0.1 we obtain a qualitatively different result. In this case 
the angular momentum transport provided by the acoustic instabil¬ 
ity is insufficiently rapid to transfer mass inward at the rate of 1 Mg 
yr -1 at which it enters the computational domain. As a result, the 
gas stagnates and collapses to form a gravitationally-unstable ring 
at a radius of 350 — 400 pc rather than ~ 100 pc as in the other 
cases. Collapse is also significantly delayed relative to the fiducial 
case, taking ~ 30 Myr to evolve to a gravitationally-unstable state 
rather than ~ 15 Myr. 

We next try varying the inflow rate M m , from 10~°' 5 — 10°' 5 
Mq yr” 1 . Figures 7 and 8 show the results. Clearly varying the 
rate at which matter enters the CMZ from larger radii in the disc 
affects the overall timescale of the evolution, but mostly in a way 
such that the behavior remains self-similar. A higher accretion rate 
simply produces a larger overall surface density at fixed time, and 
a proportionately larger mass in the gravitationally unstable region, 
but the shape of surface density and Toomre Q versus radius are 
unchanged. Similarly, in all three cases the location of the grav- 
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Figure 6. Same as Figure 4, except that we only show the mass and radius Figure 8. Same as Figure 6, but for runs with varying values of M- In as 

of the material with Q < 1 (not Q < 10), and we show results for runs indicated in the legend (in units of Mg yr — 1 . The line for M; n = 1 is the 

with varying values of ao as indicated in the legend. The line for oq = 1 same as in Figure 4. 

is the same as in Figure 4. 


i = 15Myr t—45 Myr 



Figure 7. Same as Figure 5, but with all parameters set to their fiducial 
values (Table 1) except M; n , which varies as indicated in the legend (in 
units of Mg yr — 1 ). 


itationally unstable ring is roughly the same, and rate at which 
the mass in the gravitationally-unstable region grows eventually 
asymptotes to match the rate at which matter enters the computa¬ 
tional domain. The formation of an unstable region is clearly robust 
against changes in the accretion rate. 

Our third variation is in crm, the velocity dispersion of the ma¬ 
terial entering the computational domain. This affects the rate of 
transport near the outer computational domain boundary, before 
the material advects inward far enough for its velocity dispersion 
to set by the balance between transport and dissipation. Figures 9 
and 10 show the results of simulations using <Tj n = 20, 30, 40, 
and 80 km s -1 , roughly spanning the plausible observed range for 


t = 15 Myr i =45 Myr 



Figure 9. Same as Figure 5, but with all parameters set to their fiducial 
values (Table 1) except rr; n > which varies as indicated in the legend. 

CMZ material. The results for 30, 40 and 80 km s^ 1 are clearly 
very similar qualitatively. For 20 km s -1 , the results are similar at 
times before ~ 20 Myr, but after that point the low velocity dis¬ 
persion of the material entering the computational domain leads to 
lower rates of transport at large radii, and this in turn causes mass 
to build up at large radii. Eventually the gas builds up to form a 
second gravitationally-unstable region near ~ 400 pc, in addition 
to the first one formed at ~ 100 pc. Exploration of different values 
of cri n shows that this secondary unstable region appears for values 
of (Ti n below roughly 20 — 25 km s _ . 

Finally, we vary the parameter /shape that describes how flat¬ 
tened the stellar potential is, and thus how large the gas scale height 
is at fixed gas surface density and velocity dispersion. Since the tur¬ 
bulent dissipation time is set by the scale-height crossing time, this 
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Figure 10. Same as Figure 6, but for runs with varying values of cr; n as 
indicated in the legend. The line for a = 40 km s —1 is the same as in 
Figure 4. 


t —15 Myr t —45 Myr 



100 200 300 400 100 200 300 400 

r[pc] r[pc] 


Figure 11. Same as Figure 5, but with all parameters set to their fiducial 
values (Table 1) except /shape> which varies as indicated in the legend. 

choice affects the rate at which the turbulence decays. Figure 11 
and Figure 12 show the results of this experiment. Clearly the value 
of /shape changes the precise location of the ring of gas buildup and 
gravitational instability, but does not alter the qualitative result that 
such a ring forms on timescales of ~ 10 — 20 Myr, or the default 
distribution of the gas approaching that ring. 

To summarize our findings: gas entering the CMZ tends to 
hang up in its inward flow and develop a gravitationally-unstable 
region at ~ 100 pc, where the rotation curve has a region of min¬ 
imal shear as it approach solid body. This result is very robust 
against variations in the mass accretion rate and the shape of the 
stellar potential. The former only acts as a scaling parameter, while 
the latter changes the precise location of the unstable region but 
not its general characteristics. This tendency is also robust against 


Figure 12. Same as Figure 6, but for runs with varying values of /shape as 
indicated in the legend. The line for /shape = 2 is the same as in Figure 4. 

increases in the velocity dispersion of the incoming gas and the 
rate of angular momentum transport produced by acoustic insta¬ 
bilities, and to some level of decrease in these parameters as well. 
However, if the input velocity dispersion is too low, or the rate of 
angular momentum transport by acoustic instability too small, the 
gas collapses into gravitational instability near the edge of our sim¬ 
ulation domain, either instead of or in addition to in the location of 
minimal shear at ~ 100 pc. Despite this, we can be reasonably con¬ 
fident that, over a broad range of parameters, buildup of a ~ 100 
pc unstable ring should occur on timescales of tens of Myr. 

5 DISCUSSION 

5.1 Comparison to the Structure of the Galactic CMZ 

The results presented in this paper show quantitatively how grav¬ 
itational collapse is inhibited within the ILR resonance of galac¬ 
tic centres (and of the CMZ of the Milky Way in particular). This 
is caused by efficient angular momentum transport in acoustic in¬ 
stabilities, which maintains the turbulent velocity dispersion at a 
consistently high level (a = v/3otd > 40 kms -1 ). In Kruijs¬ 
sen et al. (2014), it was proposed that the observed high turbulent 
pressure (a ~ 40 kms -1 and P/k ~ 10 s K cm -3 , see e.g.. 
Bally et al. 1988) and the correspondingly high Toomre stability 
parameter (Q ~ 10) of the CMZ gas are responsible for the fact 
that the star formation rate in the CMZ is observed to be an order 
of magnitude lower than expected from empirical scaling relations 
and star formation theory (Longmore et al. 2013a). It remained an 
open question whether the high turbulent pressure is driven by stel¬ 
lar feedback, the gas inflow along the bar, or acoustic instabilities. 
In the present paper, we find that acoustic instabilities naturally lead 
to the observed, high velocity dispersions. Our model therefore pro¬ 
vides a key missing element in explaining the currently low star 
formation rate in the CMZ. 

Next to providing an accurate description of the gas in the 
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Figure 13. Column density map of the dense (n > several 10 3 cm -3 ) gas in the Galactic CMZ as traced by NH 3 (!, 1) (Walsh et al. 2011; Purcell et at 
2012), with symbols indicating Sgr A* (green plus), star-forming clouds (red crosses; from left to right these are Sgr B2, Sgr B1, and Sgr C), and young stellar 
clusters (blue diamonds ; from left to right these are the Quintuplet and Arches clusters). The overlay provides several of our model predictions at t = 17.5 Myr 
in the plane perpendicular to the line of sight at the distance of Sgr A*. The dashed lines show the scale height profile predicted for all gas (i.e. not just the 
dense gas shown here), whereas the scale bars indicate the time necessary for the gas at each longitude to reach the gravitationally unstable region (black line), 
in units of Myr (top) as well as normalised to the local orbital time at each radius (bottom). As the gas is transported inwards, the remaining number of orbital 
times decreases more slowly than the remaining absolute time, illustrating a decrease of transport rate. As a result, the entire migration process proceeds over 
several orbital revolutions and the gas piles up just outside the gravitationally unstable region. 


CMZ on large scales, our model also reproduces the presence of 
a gravitationally unstable, ring-like stream of gas (Molinari et al. 
2011; Kruijssen, Dale & Longmore 2015). As gas loses angular 
momentum and flows towards the Galactic centre, it accumulates 
at a radius of ~ 100 pc, where the shear has a local minimum. 
This reduces the angular momentum transport rate and the driving 
of turbulence by acoustic instabilities. Consequently, the gas accu¬ 
mulates and the surface density increases, until after 10-15 Myr a 
gravitational instability develops. For the fiducial parameter set, the 
model predicts a peak surface density of E > 10 3 M & pc 2 and 
a minimum Toomre Q = 0.3-1.0 at t = 15 Myr, briefly after the 
gravitational instability has set in. These numbers are in remarkable 
agreement with the observed properties of the gas stream (cf. Ta¬ 
ble 1 of Kruijssen et al. 2014). 6 

To make a more detailed comparison, we focus on the state 
of our fiducial simulation at 17.5 Myr of evolution, which we com¬ 
pare to the observed CMZ in Figure 13; we select this time slice be¬ 
cause it has depletion times in good agreement with those observed 
in the present-day CMZ (see Section 5.2 for details). The figure il¬ 
lustrates that the gravitationally unstable region matches the radii 
where currently most of the dense gas, star-forming regions, and 
young stellar clusters in the CMZ reside (e.g., Molinari et al. 2011; 
Longmore et al. 2013b). In our model, this region covers a closed 
circular ring by definition, because we assume axisymmetry on the 
sub-kpc scales under consideration. However, the gas itself is ex¬ 
pected to follow a possibly eccentric, stream-like structure within 
this unstable radius interval, which must be open-ended due to the 


6 In our model, the typical values of E and Q within the gravitation¬ 
ally unstable radius interval are primarily set by M and cr; n (see Sec¬ 
tion 4.2). For good agreement with the observations, the model requires 
M ~ 1 Mg yr —1 and o-; n > 30 km s —1 . The former is consistent with 
other theoretical constraints (e.g., Crocker 2012), whereas the latter of these 
is consistent with the observed velocity dispersions mentioned above. 


extended nature of the stellar mass distribution. The radial range 
across which the gravitational instability takes place in our param¬ 
eter survey is 60 < R/ pc < 120, which indeed matches the radii 
covered by the eccentric orbital model for the CMZ gas stream by 
Kruijssen, Dale & Longmore (2015). 

At radii beyond the gravitationally unstable region (i.e. 150 < 
R/ pc < 500 or 1° < |(| < 3.5° in projection), dense clouds 
may exist, but in our model they have high velocity dispersions and 
are supervirial due to acoustic instabilities (a v ir ~ P tur b / -Pgrav = 
a 2 /nGHgT, = 10-100). Therefore, their global collapse is inhib¬ 
ited and the little star formation activity they are expected to have 
will result from stochastic turbulent motion (cf. Bania’s Clump, 
Bally et al. 2010). For the CMZ model presented here, star forma¬ 
tion theories (e.g., Krumholz & McKee 2005; Padoan, Haugbplle & 
Nordlund 2012; Federrath & Klessen 2012; Hennebelle & Chabrier 
2013) predict low star formation efficiencies per free-fall time of 
eff ~ 0.001 outside the gravitationally unstable region. This is par¬ 
ticularly true if the turbulence is primarily solenoidal (Federrath & 
Klessen 2012), as might be expected since it is ultimately driven by 
shear. 

As illustrated in Figure 13, the high velocity dispersions re¬ 
sult in an increase of the gas scale height with radius. Depending 
on their Galactic longitudes, it takes these highly-turbulent clouds 
anywhere between 2-15 Myr to migrate inwards to the gravitation¬ 
ally unstable region at R ~ 100 pc (cf. the third panel of Figure 3, 
which shows inward radial velocities of — v v = 20-30 km s _ ), 
corresponding to ~ 2 full orbital revolutions. This estimate as¬ 
sumes that the clouds all reside at the same distance - for non-zero 
displacements along the line of sight, the migration time-scales are 
even longer. Direct-infall models for the gas inflow onto the CMZ 
(e.g., those in which the *i and X 2 orbits play a major role) predict 
much shorter migration time-scales of 0.5-5 Myr (e.g. Sofue 1995; 
Bally et al. 2010). Comparing observed cloud properties as a func¬ 
tion of longitude with the expected evolutionary time-scales may 
therefore be a way of discriminating between direct-infall models 
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Table 3. Predicted dense gas fractions in the Galactic CMZ. 


Longitude 

./dense,4 

(> 10 4 cm -3 ) 

/dense, 5 

(> 10 5 cm -3 ) 

/dense,6 
(> 10 6 cm -3 ) 

1 1\ = 2-3° 

0.03 

0.003 

0.0002 

1*1 ~ 1-5° 

0.05 

0.006 

0.0005 

1*1 ~1° 

0.08 

0.015 

0.0015 

|t| < 0.7° 

0 .3-1.0 

0.04-0.2 

0.005-0.01 


These values apply to our fiducial model (Figure 3). While the overall 
trend of increasing dense gas fractions towards lower longitudes 
is a robust prediction of the model, the normalisation of that trend 
will depend on the adopted parameters (cf. Section 4.2). As these 
parameters become better constrained by new observations, so will the 
dense gas fractions provided here. 

and the dynamical evolution model presented here. Other tests of 
the model can be performed by comparing the radial profiles of 
e.g., the gas surface density and velocity dispersion from Figure 3 
to observed maps of the CMZ. 

As a result of the inward transport of the gas and the corre¬ 
sponding increase of the gas pressure towards lower Galactic lon¬ 
gitudes, our model predicts that the dense gas fraction increases 
strongly towards the Galactic Centre. To quantify this prediction, 
we consider gas near the disc midplane, where the mean density is 
p = E/2 H g . We then assume that the gas has a log-normal volume 
density PDF as expected for supersonically turbulent, isothermal 
media (e.g. Vazquez-Semadeni 1994; Padoan, Nordlund & Jones 
1997; Krumholz & McKee 2005), and calculate the mass fraction 
above a certain minimum density. 7 The resulting dense gas frac¬ 
tions are listed at four different Galactic longitudes in Table 3. 
Interestingly, the dense gas fractions change only slowly in the 
outer CMZ, where acoustic instabilities inhibit gravitational col¬ 
lapse. However, they increase rapidly as the gas becomes gravita¬ 
tionally unstable at |2| < 1°. Designating all gas with densities 
n > 10 4 cm -3 as ‘dense’, we find that the dense gas fraction in¬ 
creases from a few per cent in the outer CMZ to nearly unity in the 
100-pc ring (as was also found by Longmore et al. 2013a). These 
fractions decrease strongly as the minimum density used to define 
‘dense’ gas increases. Interferometric observations of high-critical 
density gas tracers using ALMA or the ongoing SMA Legacy Sur¬ 
vey of the CMZ (Keto et al. in preparation; Battersby et al. in prepa¬ 
ration) are ideally suited to test these predictions. 

5.2 Time-evolution of the Galactic CMZ 

It is important to reiterate that our model predicts an episodic cy¬ 
cle for the gas content and star formation activity of galactic cen¬ 
tres. The development of a gravitationally unstable region marks 
the beginning of the end: once the gas collapses and forms stars, 
the residual gas is expelled by feedback. The current gas mass of 
the 100-pc stream in the CMZ is M ~ 10‘ M 0 , with a density 
of n > 10 4 cm -3 and a correspondingly short free-fall time of 
tg < 0.34 Myr, which is much shorter than the time required to 
grow the current gas reservoir (t acc ~ M/M[ n ~ 10 Myr for our 

' The density PDF can be altered somewhat depending on the ratio of 
solenoidal to compressive modes in the driving force (Federrath, Klessen 
& Schmidt 2008), the strength of magnetic fields (Molina et al. 2012), and 
deviations from isothermality (Federrath & Banerjee 2015), but for sim¬ 
plicity we perform this calculation using the results for an isothermal, non¬ 
magnetic medium with mixed solenoidal-compressive driving. 


fiducial model). As a result, the gravitationally unstable gas reser¬ 
voir is expected to be depleted by a combination of rapid star for¬ 
mation and feedback, after which the accretion cycle described by 
our model starts from the beginning. 

Given the episodic nature of the system's evolution, a more 
detailed comparison between model and observations requires us 
to assess at which point along the cycle the CMZ currently re¬ 
sides. The present star formation rate in the CMZ is only ~ 
0.05 M 0 yr” 1 (Longmore et al. 2013a; Koepferl et al. 2015) and 
the mass outflow rate is estimated at ~ 0.5 M 0 yr” 1 (correspond¬ 
ing to a mass loading factor of i ~ 10, see e.g., Crocker 2012). 
The sum of both (M out ~ 0.6 M 0 yr -1 ) is lower than our fidu¬ 
cial mass inflow rate of Mj n ~ 1 M 0 yr” 1 , and also lower than 
the independent estimate by Crocker (2012), who finds 2 a limits of 
0.4 < Mi n /M 0 yr” 1 < 1.8. It is therefore most probable that the 
gas mass in the CMZ is presently increasing. This places the CMZ 
at the moment prior to the starburst, but after the gravitational in¬ 
stability has started to develop. 

Comparing the current mass of the 100 pc gas stream (M ~ 
10' M 0 ) to the time-evolution of the unstable mass in our fiducial 
model (see Figure 4), the CMZ appears to reside some t ~ 20 Myr 
after its last major starburst. However, this is an upper limit, as it 
assumes that (1) the entire gas reservoir needed to be regrown out 
to R > 400 pc, and (2) that all of the mass in the gas stream is al¬ 
ready gravitationally unstable. In the more likely scenario that the 
previous starburst did not clear out any of the gas beyond the un¬ 
stable region, the evolution towards gravitational instability during 
the first ~ 10 Myr is skipped, because accretion can resume with¬ 
out delay. The growth phase of the gravitational instability then re¬ 
mains, resulting in a time interval of ~ 10 Myr to accumulate the 
present-day gas mass of the 100 pc stream. If the previous starburst 
did not clear out all of the gas in the unstable region, this time in¬ 
terval will be further reduced. We therefore conclude that the CMZ 
is currently 5-10 Myr post-starburst. 8 

The obvious next question to ask is when the CMZ will en¬ 
ter the next starburst phase. To address this point, we consider 
our model at t = 17.5 Myr, shortly after the initial formation 
of the gravitationally-unstable mass reservoir. If the properties of 
these gravitational instabilities match those observed in the CMZ, 
their growth time-scales provide insight in the CMZ’s immediate 
evolution. Figure 14 shows the length-scale of the acoustic and 
gravitational instabilities as well as the growth time-scale of the 
fastest-growing mode as a function of radius in the CMZ, at times 
t = {10,17.5} Myr. For the acoustic instabilities, the figure gen¬ 
eralises the initial result from Figure 2 to the evolved CMZ model, 
showing that the growth time-scale of these instabilities is typically 
less than an orbital time and that they develop rapidly enough to 
drive angular momentum transport and turbulence. 9 


8 The main observational relics of this previous starburst are likely the 8- 
kpc Fermi bubbles that extend from the CMZ (Su, Slatyer & Finkbeiner 
2010). In combination with their spatial extent, the time since the starburst 
given here constrains the outflow velocity of the bubbles to be « 0 ut = 800 
1600 km s” 1 , which is remarkably consistent with previous estimates from 
geometric considerations (n ou t ~ 1000 km s”, Carretti et al. 2013). 

9 The characteristic length-scale of the instabilities is short, at A < 
0.01 pc. However, it is important to recall that this wavelength charac¬ 
terizes the radial size in the tight-winding approximation; that is, it is the 
characteristic width of the m = 2 spiral features, and tells us nothing about 
the pitch angle of the spiral, other than the fact that it must be small since 
the Montenegro, Yuan & Elmegreen (1999) dispersion relation is derived in 
the limit where it is. Moreover, we emphasize that the fastest growing linear 
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Figure 14. Unstable wavelengths for the acoustic (shaded area) and grav¬ 
itational (hashed area) instabilities in the fiducial model, with the fastest- 
growing mode colour-coded by the growth time-scale in units of the orbital 
time (see colour bar). The top panel shows the radial profile prior to the on¬ 
set of the gravitational instability at t = 10 Myr, whereas the bottom panel 
shows the “current-day” state at f = 17.5 Myr. 


The development of the gravitational instabilities around R ~ 
100 pc differs from the acoustic instabilities. Figure 14 demon¬ 
strates that the growth time-scale is initially longer than an orbital 
time. 10 As the instability develops, the growth time-scale decreases 
substantially to f gro wth < forb, and two modes with different wave¬ 
lengths dominate: 

(i) At radii R = 80-100 pc, the instability length-scale is 
A ~ 200 pc. This long-wavelength instability is expected to drive 
large-scale asymmetries in the distribution of the gas in the CMZ, 
as it covers a significant fraction of a circular orbit at these radii. 
This matches the well-known observation that most of the gas in 
the CMZ resides at positive Galactic longitudes (see Figure 13 and 
Bally etal. 2010). 

(ii) At radii R = 70-80 pc, the fastest-growing mode is a fac¬ 
tor of several shorter, with A ~ 30 pc. This short-wavelength in¬ 
stability matches the estimated Jeans length in the 100-pc stream 


mode is not necessarily the dominant one in the non-linear regime, and thus 
the dominant mode that we would expect to observe. An obvious example 
is Rayleigh-Taylor instability, for which the growth rate in the linear regime 
scales as A 1 ' 2 , so that small-scale modes are fastest, but both simulations 
and experiments show that long-wavelengths mode dominate in the non¬ 
linear regime (e.g., Dimonte et al. 2004). In the absence of full simulations 
of the acoustic instability, we cannot reach any strong conclusions about 
which modes will dominate in real galaxies. 

10 This supports the scenario put forward by Longmore et al. (2013b) and 
Kruijssen, Dale & Longmore (2015) in which the collapse of clouds on the 
100 -pc stream is triggered by a tidal compression during pericentre passage. 
Because the gas on the stream initially requires several orbital revolutions 
to become gravitationally unstable, the recurrent tidal perturbations at peri¬ 
centre are statistically likely to give the final nudge into collapse as the 
turbulent energy is gradually dissipating. 


(~ 30 pc, Kruijssen, Dale & Longmore 2015), as well as the ob¬ 
served typical separation length of the individual clouds condens¬ 
ing out of the stream and the wavelength of line-of-sight velocity 
oscillations (Henshaw et al. in preparation). 

Given how well these instabilities match the observed density fluc¬ 
tuations in the CMZ, their growth time-scales may be taken as an 
indication of how rapidly the CMZ clouds can start their collapse 
towards the next starburst phase. 

Both of the above gravitationally unstable modes have 
t growth < forb ~ 3.7 Myr. Therefore, the extant condensations 
should start collapsing within the next few Myr, and the precise 
time delay till the next starburst within the gravitationally unstable 
region will depend on the number of free-fall times needed to fun¬ 
damentally change the ratio between gas mass and star formation 
rate. We use the observational estimate that star formation efficien¬ 
cies of e s f ~ 0.1 are reached before a star formation event dis¬ 
perses the residual gas (Lada & Lada 2003; Federrath & Klessen 
2013). The associated time-scale is t s { = tge s f/eg, where eg is 
the star formation efficiency per free-fall time. In the Galactic disc, 
eg ~ 0.01 (Krumholz & Tan 2007; Krumholz, Dekel & McKee 
2012; Federrath 2013), but it is unclear if this holds in the CMZ 
as well. 11 Using orbital modelling, Kruijssen, Dale & Longmore 
(2015) estimate that the evolutionary time difference between the 
'Brick’ (little star formation) to Sgr B2 (high star formation activ¬ 
ity) is Af ~ tg ~ 0.4 Myr, whereas the star formation efficiency 
in Sgr B2 is of the order a per cent (Bally et al. 2010, assuming a 
few 100 Mq per ultracompact Hll region). As a result, eg ~ 0.01 
in collapsing CMZ clouds too. Combining the above numbers, we 
estimate that the next starburst in the gravitationally unstable region 
will be reached in ~ 5 Myr or 1-2 orbital revolutions. 


5.3 A Cartoon Model with Star Formation and Feedback 


We are now in a position to develop a toy model for the overall 
behaviour of star formation in central molecular zones, both those 
of the Milky Way and of similar galaxies, and how they evolve in 
the observational plane of star formation rate surface density ver¬ 
sus gas surface density. Our goal here is not to precisely reproduce 
the properties of the Milky Way’s CMZ. Instead, it is to develop a 
simple cartoon picture for how regions such as the CMZ evolve in 
time. To do so, we must extend our model with a simple treatment 
of star formation and feedback. 

To this end, we modify our fiducial model by adding a term 
—S* on the right-hand side of equation 2. We compute the star 
formation rate as 

St = fcieg — (18) 

tg 

where f c i is the fraction of gas in dense clouds, eg is the star forma¬ 
tion rate per free-fall time in those clouds, tg is the free-fall time. 
We compute these quantities as f c i = 1/2, 


eg = 0.01e"“ vir 

P/H g 

(tt/2)GS 2 


(19) 

( 20 ) 


11 When globally-averaged surface densities are adopted, the CMZ is con¬ 
sistent with galactic star formation relations and models assuming eg ~ 
0.01 (Yusef-Zadeh et al. 2009; Longmore et al. 2013a; Kruijssen et al. 2014; 
Salim, Federrath & Kewley 2015), but this approach ignores persistent (and 
hence physical) substructure. Therefore, it does not directly constrain the 
true value of eg in the CMZ. 
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te = yf&fp 

E 

P ~ 2If,' 


( 21 ) 

( 22 ) 


Physically, these expressions amount to saying that half the mass 
at any given radius is in clouds as opposed to in a diffuse medium, 
roughly consistent with observations of circumnuclear starbursts 
(e.g., Rosolowsky & Blitz 2005). Within those clouds, star for¬ 
mation proceeds at the observed rate of ~ 1% of the mass per 
free-fall time in virialized gas, again consistent with observations 
(Krumholz & Tan 2007; Krumholz, Dekel & McKee 2012; Fed- 
errath 2013; Salim, Federrath & Kewley 2015), but that the star 
formation rate drops off exponentially in super-virial gas (that with 
Ovir > 1). Our definition of the virial ratio is such that a v i r —> 1 
when gas self-gravity dominates over stellar gravity in equation 6, 
indicating that the gas is self-gravitating. 

When we re-run the fiducial case with this added term, the 
behaviour is qualitatively unchanged, except that rather than gas 
accumulating indefinitely in the unstable region, the ring instead 
reaches a steady state where rate of star formation in the ring bal¬ 
ances the rate of gas transport into it. This behaviour is quite insen¬ 
sitive to the exact star formation recipe we adopt. In reality, such a 
steady state is unlikely to be stable. The star formation rate surface 
density in the steady state is so high that the radiation flux reaches 
~ 10% of the Eddington luminosity, at which point radiation pres¬ 
sure alone should drive rapid mass loss, even without the aid of 
supernovae (Thompson & Krumholz 2014). To construct our toy 
model, we therefore add an artificial truncation to the gas accumu¬ 
lation and star formation. After 22.5 Myr of evolution (chosen to be 
5 Myr after the ‘current’ CMZ as discussed in Section 5.2), when 
the star formation rate has reached ~ 0.1 M 0 yr -1 , we simply 
reduce the star formation rate linearly to zero, and the gas surface 
density linearly to 10% of its previous value, over a time of 4 Myr, 
roughly the lifetime of a massive star, in every computational cell 
for which Q < 3. 12 We note that this timescale is selected to give 
a qualitative picture, and is not quantitatively calculated. 

With this quite crude cartoon model of feedback and gas ex¬ 
pulsion, the overall evolution of the star formation rate, gas mass, 
and depletion time is as shown in Figure 15. In this plot we show 
values computed both using all the material in the inner 250 pc, 
and values computed only for the material within the ring from 
80 — 90 pc ring where the gas surface density and star forma¬ 
tion peak. Our motivation for separating these two cases is that 
the former is roughly analogous to what would be observed in a 
survey with ~ 0.5 kpc resolution (e.g., HERACLES, Leroy et al. 
2013), targeting a galaxy > 1 Mpc away using PdBI or a similar 
instrument; in our own Galaxy, it is roughly equivalent to consid¬ 
ering all the material with Galactic longitude |2| < 3°. The latter 
is roughly analogous to considering the material within |2| < 1° 
in the Milky Way, or to the resolution that is possible in external 
galaxies using ALMA. We also construct a realistic “observed” star 
formation rate by convolving the true star formation rate with the 
time-dependent ionizing luminosity of a stellar popula¬ 

tion that fully samples the IMF. This correction is important be¬ 
cause galactic centre star formation rates are usually measured with 
Ha or similar ionization-based star formation rate indicators, and 
these effectively average the star formation rate over a time that is 


12 To avoid discreteness noise, we use the same trick mentioned above, 
whereby we interpolate Q and all other quantities onto a higher resolution 
grid and apply this procedure to the interpolated data. 


0 

0 

„ 0 
k 0 

sf 0 

e£ 0 

UL 

ui Q 
0 
0 


.40 

35 

.30 

25 

.20 

15 

10 

.05 

.00 

20 

"0 

a 15 

7 10 

a 

5 

0 

10 1 


>. 10 u 


10 ” 


10 ” 



10 


15 

t [Myr] 


20 


25 30 


Figure 15. Star formation rate (top panel), total gas mass (middle panel), 
and gas depletion time (/.,] ep = M gas /SFR, bottom panel) for our fiducial 
run with star formation added. In all panels, quantities indicated by thin 
lines show true, instantaneous values, while those indicated by thick lines 
show values that would be inferred from observations using an ionization- 
sensitive star formation rate indicator such as Ha. Dashed lines show a run 
with no gas removal, while solid lines show the results where we use our 
crude gas expulsion model, whereby the star formation rate declines linearly 
to zero over a 4 Myr time period starting at 22.5 Myr (see main text for 
details). In the middle and bottom panels, the blue lines indicate values over 
the inner 250 pc, as would be observed with ~ 0.5 kpc resolution, while 
green lines indicate values measured only in the ring between 80 and 90 pc 
from the Galactic Centre, as could be observed with ~ 10 pc resolution. 
We do not include the green lines for 10 pc-resolution in the upper panel 
because they are essentially identical to the blue lines, since nearly all of 
the star formation takes place within the ring. 


non-negligible compared to timescales in the model. We compute 
Q(H°,f) using the slug code (da Silva, Fumagalli & Krumholz 
2012; Krumholz et al. 2015) using Geneva non-rotating evolution¬ 
ary tracks at Solar metallicity (Ekstrom et al. 2012). 

What does this evolutionary cycle look like when plotted on 
the usual star formation relations, which relate the star formation 
rate per unit area to the gas surface density, or surface density nor¬ 
malized by orbital time? To answer this question, we again con¬ 
sider observations with both ~ 0.5 kpc and ~ 10 pc resolution. 
For the former, we take the area within a radius of 0.25 kpc (i.e. 0.2 
kpc 2 ) and consider all the material within this radius (cf. Kruijs¬ 
sen et al. 2014), while for the latter we use the area of an annulus 
from 80 — 90 pc, which is where essentially all the star formation 
is located - this corresponds to roughly the central 1° in Galactic 
longitude. For the orbital time, we use the orbital time for our fidu¬ 
cial rotation curve evaluated at the outer edge of the disc for the 
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~ 0.5 kpc-resolution case, and the orbital time evaluated at 90 pc 
for the ~ 10 pc-resolution case. 

We must also choose when to start and stop the clock in or¬ 
der to make this comparison. As noted above, our initial condition 
contains an unrealistically small gas mass, since even a strong star- 
burst in the unstable region seems unlikely to expel material that is 
~ 400 pc from a galactic centre. With this consideration in mind, 
we choose our zero of time to be 14 Myr after the start of the fidu¬ 
cial run, which marks the moment when a region with Q = 1 first 
appears. With this choice, and using our crude model to represent 
star formation feedback, gas accumulates for 8.5 Myr with a slowly 
increasing star formation rate. After that point gas expulsion begins 
and star formation slows down. The expulsion process lasts until 
12.5 Myr, and thereafter the stellar population fades and the ob¬ 
served star formation rate (which lags the true one) declines, until 
gas begins re-accumulating and the cycle repeats. The full cycle in 
this cartoon model requires 15-20 Myr. 

Figure 16 shows the cycle as it appears in the star forma¬ 
tion relations. Clearly at either spatial resolution the variation of 
the CMZ’s gas depletion time is dominated by the fluctuating star 
formation rate. The gas mass within the entire CMZ is constant to 
within a factor of 2 — 3. The mass inside the unstable ring varies 
more widely, but the amount of variation is not well-determined, 
since it depends strongly on our assumptions about how much mass 
is expelled when the ring begins a starburst. Despite this uncer¬ 
tainty, it seems likely that the gas depletion time of the gravita¬ 
tionally unstable, inner CMZ will vary by at least two orders of 
magnitude, while the depletion time of the CMZ as a whole will 
show ~ 1 dex variations. This difference arises because on large 
scales, the independent, outer gas reservoir is included in the de¬ 
pletion time measurement, which stabilises the small-scale fluctua¬ 
tions (Kruijssen & Longmore 2014). 

Focussing on the left panel of Figure 16, and comparing the 
behaviour in our toy model to the star formation relations observed 
for entire galaxies (various lines in Figure 16), we see that a CMZ 
measured on ~kpc scales spends roughly half its time with a star 
formation rate below what would be expected based on its gas sur¬ 
face density, and about half the time at a higher star formation 
rate. The inner CMZ (i.e., the 100-pc stream) might seem to spend 
more time above the mean star formation relation, but this is some¬ 
what misleading. The Bigiel et al. (2008) relation, marked by the 
red dashed line in Figure 16, represents a fit to disc galaxies with 
E < 10 2 Mq pc -2 , and does not describe the elevated star for¬ 
mation rates of galaxies with higher surface densities. Compared to 
the star formation relations observed at these higher surface densi¬ 
ties, a CMZ observed at 10 pc-resolution also spends roughly half 
its time below the mean star formation rate, and half its time above. 
For ~ 0.5 kpc-resolution observations, the scatter about the con¬ 
ventional star formation relations is ~ 1 — 1.5 dex, while at ~ 10 
pc resolution it is as much as ~ 2 — 2.5 dex, though the exact value 
depends on the details of gas expulsion. The present-day CMZ of 
the Milky Way is in the phase of its evolution when it lies below 
the mean star formation relations. 

The above comparison changes when dividing the gas surface 
density by the orbital time (right-hand panel in Figure 16). In that 
case, the extremely short orbital time of a CMZ places it below all 
star formation relations except during its starbursts. This holds for 
both the inner CMZ and the entire CMZ. The division by the or¬ 
bital time therefore provides an efficient way of highlighting galac¬ 
tic centres at a star formation minimum, irrespective of the spatial 
resolution. 

Finally, we note that in all four discussed cases (two spatial 


scales in two panels of Figure 16), the cycle-averaged CMZ agrees 
with Daddi et al. (2010)’s observed star formation relations for star- 
burst galaxies to within a factor of 3. Without any significant evo¬ 
lution of the rotation curve or the gas inflow rate (either by secular 
processes or due to external perturbations), the CMZ should con¬ 
tinue to evolve through the cycle shown in Figure 16. If the condi¬ 
tions do change fundamentally, then this will affect the time-scales 
of the different phases. The cycle is not expected to shut off alto¬ 
gether, because there should always be a radius interval where the 
rotation curve has a near-solid body part and gas accumulates. 

5.4 Predictions for the Centres of External Galaxies 

The model presented here provides the first self-consistent theory 
that quantitatively reproduces the main features of the CMZ of the 
Milky Way. It is therefore desirable to make predictions that hold 
for galactic centres of barred spiral galaxies in general. 

In summary, our model consists of the following critical in¬ 
gredients. 

(i) The galaxy has an ILR to which gas is supplied by the bar 
at sufficiently high velocity dispersions (ctid > 10-15 km s -1 ) to 
drive angular momentum transport and thus prevent the gas inflow 
from stalling. 

(ii) For most of the radial range within the ILR, the rotation 
curve is near-flat so that the shear is high (1 — /3 > 0.8). 

(iii) As a result of points (i) and (ii), acoustic instabilities de¬ 
velop in the gas after it enters the ILR. These instabilities drive 
efficient angular momentum transport, leading to inflowing gas 
with high turbulent pressures (P/k > 10 5 K cm -3 ), extreme 
gravitational stability (Q > 10), and a low star formation rate 
(idep ~ 10 Gyr). 

(iv) Within the ILR, there is a radius where the rotation curve 
transitions from near-flat to near-solid body, causing a minimum 
in the shear (1 — /3 < 0.4). The angular momentum transport and 
mass inflow both stall at this radius. 

(v) At this radius, the inwards mass flux is sufficient for the ac¬ 
cumulating gas to eventually exceed the stellar density (p > p*). 

(vi) As a result of points (iv) and (v), gravitational instabilities 
develop in the growing gas reservoir near the shear minimum. At 
the high gas densities required by point (v), these instabilities drive 
the gas into rapid free-fall (tfj < 1 Myr), resulting in a starburst 
(fdep <0.1 Gyr). 

These conditions are commonly satisfied in the centres of barred 
spiral galaxies (e.g. Elmegreen, Elmegreen & Eberwein 2002; Jo- 
gee et al. 2002; Martini et al. 2003; Jogee, Scoville & Kenney 2005; 
Leroy et al. 2008; Sandstrom et al. 2010; Nesvadba et al. 2011 ; Sani 
et al. 2012). As a result, the physical processes described in this 
paper and the qualitative cycle of Figure 16 (also see Figure 6 of 
Kruijssen et al. 2014) should be common in other galactic centres. 

While our model may apply qualitatively to other galaxies, it 
is unclear if the quantitative predictions hold. Given the close re¬ 
lation between the Galactic rotation curve and the predictions of 
our model, it is to be expected that the details of extragalactic ro¬ 
tation curves lead to quantitatively different radial profiles of, e.g., 
the gas surface density, velocity dispersion, and mass inflow rate. 
However, the duty cycle from quiescence to starburst activity and 
back may be relatively unaffected. The ILR of more massive galax¬ 
ies may reside at larger radii (which at first sight implies longer gas 
accretion time-scales), but the deeper gravitational potential and the 
correspondingly increased gas pressures will accelerate the angular 
momentum transport by acoustic instabilities (which decreases the 
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Figure 16. Evolution of a CMZ in the planes of star formation surface density Ssfr versus gas surface density S ( left panel) and EgpR versus gas surface 
density divided by orbital time E/l or b. Blue lines show the system as it would be observed with ~ 0.5 kpc resolution, while green lines show the system at 
~ 10 pc resolution (see text for details). Circles mark various critical points in the evolution: the start of gas accumulation at 0 Myr (corresponding to 14 Myr 
of evolution in the fiducial run), the onset of gas expulsion at 8.5 Myr, the end of gas expulsion at 12.5 Myr, and complete fading of the ionizing flux from 
the stellar population at ~ 17 Myr. The stars, labelled “CMZ”, mark the rough point in this evolutionary cycle at which the Milky Way’s CMZ resides, and 
corresponds to 3.5 Myr after time 0 in the cycle. Diamonds indicate the time-average over the full cycle. For comparison, we also show the star formation 
relations of Kennicutt (1998, black dotted line), Bigiel et al. (2008, red dashed line), and Daddi et al. (2010, magenta dot-dashed line). For Kennicutt (1998), 
we have multiplied the star formation rate in the original fit by 1.6 to adjust from Kennicutt’s adopted Salpeter (1955) IMF to a Chabrier (2005) one. For Daddi 
et al. (2010), the two lines shown in the left panel correspond to Daddi et al.’s separate fit to discs ( lower line) and starbursts (upper line). 


gas accretion time-scale again). As a result, it is possible that the 
duty cycle predicted for the CMZ (see Figure 16) has a broader 
application than perhaps naively expected. From an empirical per¬ 
spective, it is therefore interesting to make a simple comparison to 
the population statistics of other galactic centres. 

A key result of the discussion in Section 5.3 is that the details 
of the star formation cycle and its position relative to empirical star 
formation relations depend on the spatial resolution (compare the 
blue and green tetragons in Figure 16). Low-resolution observa¬ 
tions should find a smaller variation of the gas depletion time than 
high-resolution observations. Leroy et al. (2013, Figure 13) find a 
~ 1 dex scatter of the gas depletion time in the central R < 0.5 kpc 
of 30 nearby disc galaxies from the F1ERACLES survey. A simi¬ 
lar range is found for other galaxy samples (e.g., Sakamoto et al. 
1999; Jogee, Scoville & Kenney 2005; Hsieh et al. 2011; Sani et al. 
2012; Saintonge et al. 2012; Fisher et al. 2013). In addition, the 
galactic centres in the HERACLES sample exhibit a ~ 0.3 dex de¬ 
crease of the gas depletion time relative to galactic discs (Leroy 
et al. 2013). Both the observed scatter and deviation of the deple¬ 
tion time present a remarkably good match to our predictions for 
the same area (|Z| < 3° and 0.5 kpc resolution). 

In view of the good agreement between our model and the 
galactic centres in the HERACLES sample, we can make a predic¬ 
tion for the relative frequency of quiescent, normal and starbursting 
galactic centres across the population of barred spiral galaxies. Of 
course, this assumes that our CMZ model describes the duty cy¬ 
cle despite obvious differences in rotation curves and other galaxy 


properties. In turn, large observational samples of galactic centres 
can be used to test this assumption. 

We estimate the relative occurrence rates by using the time- 
scales associated to each of the phases in the evolutionary cycle 
from Section 5.3. For 0.5 kpc-resolution observations, we predict 
that ~ 1/3 of all galaxies have a ‘normal’ gas depletion time-scale 
of ~ 2 Gyr (the Bigiel et al. (2008) line in Figure 16) within a 
factor of 2, whereas ~ 2/3 have elevated star formation activity 
and shorter gas depletion times; a small fraction have depressed 
star formation rates. This is again consistent with the distribution 
of gas depletion times found by Leroy et al. (2013). By contrast, 
for high-resolution (~ 10 pc) observations, we predict that the de¬ 
pletion times will almost always be shorter than the ~ 2 Gyr found 
at larger galactic radii, and that ~ 1/2 of the time they will also 
lie above the star formation relations of Kennicutt (1998) or Daddi 
et al. (2010) for “normal” (i.e., non-starbursting) galaxies. This is 
a fundamentally different distribution than at low resolution, and 
shows that upcoming ALMA observations provide a unique oppor¬ 
tunity to test our predictions. 

We reiterate one of the conclusions from Section 5.3, that di¬ 
viding the gas surface density by the orbital time (which is natu¬ 
rally short in galactic centres) shifts galactic centres below the re¬ 
lation followed by galaxy discs at any spatial resolution at almost 
all times (see the right-hand panel of Figure 16). Only in outburst 
do these star formation relations approach the usual ones found at 
larger galactic radii. This provides a way of highlighting galactic 
centres with low star formation rates. 

More direct tests of our model will be enabled by high- 
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resolution ALMA observations of nearby galactic centres. Such 
observations provide the necessary input for our model (i.e., high- 
resolution rotation curves) as well as the predicted observables 
(e.g., radial profiles of the gas surface density, velocity dispersion, 
mass inflow rate, dense gas fraction, etc.). Next to building up a 
statistically representative sample of galactic centres, such observa¬ 
tions have the additional advantage that the centres of other barred 
spiral galaxies may not be at the same phase of the evolutionary 
cycle as the Galactic CMZ. That way, they provide the best oppor¬ 
tunity to characterise the relevant physics during all evolutionary 
phases and to derive an accurate duty cycle for star formation ac¬ 
tivity in the centres of barred spiral galaxies. 


5.5 Relation to Previous and Future Theoretical Work 

We end this discussion by commenting on how the model we 
present here relates to other theoretical work on the dynamics of 
gas in galactic CMZs. There have been a number of previous hydro- 
dynamic simulations of such regions (e.g., Rodriguez-Fernandez & 
Combes 2008; Baba, Saitoh & Wada 2010; Pettitt et al. 2014; Cole 
et al. 2014; Sormani, Binney & Magorrian 2015). Most of this work 
has been focused either on explaining the observed kinematics of 
gas near the Galactic center, or on studying how a bar can drive 
mass inward and create circumnuclear structures. The picture that 
these simulations present is that a bar, by creating self-intersecting 
Xi orbits, leads the gas to shock and plunge inward, whereupon it 
finds itself on X 2 orbits at distances of a few hundred pc from the 
Galactic center. This behavior is roughly consistent with the pre¬ 
dictions of analytic models (Binney et al. 1991). 

With the exception of Sormani, Binney & Magorrian (2015), 
who reach 5 pc resolution, none of the previous works of which we 
are aware have enough resolution to model the ~ 100 pc size-scale 
structures on which we focus here. Because they are focused on 
larger scales, then tend to have resolutions of ~ 50 pc at best. 13 
Sormani, Binney & Magorrian find that the gas on x\ orbits that 
shocks does subsequently settle into a disc-like structure a few hun¬ 
dred pc in size, with the exact settling radius in their simulations 
depending on their assumed gas equation of state. Given the rela¬ 
tively small wavelength of the acoustic instability (cf. Figure 14), 
they unfortunately cannot resolve it, but in principle similar sim¬ 
ulations at higher resolution should be able to model the scenario 
that we propose here. 

More broadly, future simulations would provide a valuable test 
of the central assumption in our picture. We have assumed that, 
once gas is transferred inward to X 2 orbits by shocks, we can ap¬ 
proximate its behavior as an axisymmetric disc subject to pertur¬ 
bations from the bar. This approach emphasizes the gaseous nature 
of the flow in the inner parts of a galaxy. In contrast, many previ¬ 
ous authors have instead chosen to visualize the gas in this region 
as consisting of effectively-collisionless clouds moving on ballistic 
orbits (e.g., Binney et al. 1991; Sofue 1995), which plunge into the 
central black hole on an orbital timescale. If this latter picture is 
closer to the truth, then our disc-based picture is called into ques¬ 
tion. While Sormani, Binney & Magorrian (2015)’s simulation cer¬ 
tainly does suggest the formation of a disc-like structure, the ques- 


13 Rodriguez-Fernandez & Combes (2008) use a sticky-particle method 
with a ~ 5 pc collision radius, but the point still holds - because they 
are not actually simulating hydrodynamics, their method is not capable of 
representing an accretion disc. 


tion of how well gas inside the x\ region can be approximated as a 
perturbed disc remains unsettled. 


6 CONCLUSION 

In this paper we present the first global model for gas transport and 
star formation in the centres of barred spiral galaxies, focusing on 
the Milky Way’s Central Molecular Zone (CMZ) as a paradigmatic 
example. The model is based on a few simple ingredients. First, 
the rotation curves of galaxies, which are generally close to flat at 
large radii, must eventually turn over and approach solid body near 
galactic centres, producing a minimum in the local shear. In the 
CMZ, this minimum occurs ~ 100 pc from the Galactic centre. 
Second, gas that orbits in the flat rotation curve part of the poten¬ 
tial will be unstable due to acoustic instabilities pumped within the 
bar’s inner Linbdlad resonance (ILR). These instabilities will trans¬ 
port angular momentum and cause mass to flow inward, until the 
flow stagnates in the low-shear region where angular momentum 
transport is suppressed. In our axisymmetric models this produces 
a ring-like structure, though in a real galaxy the ring is likely to be 
partially rather than completely filled. Third, within the region of 
rapid infow, the velocity dispersion driven by the instability render 
the gas highly-supervirial, and thus unable to undergo collapse and 
star formation. In contrast, the accumulating gas in the stagnation 
region must eventually become self-gravitating and undergo vigor¬ 
ous star formation, likely leading to a blowout that removes much 
of the gas. The cycle will then begin again as new mass accumu¬ 
lates. 

The scenario we propose, which is supported by numerical 
calculations of disc evolution performed with the VADER code 
(Krumholz & Forbes 2015), naturally explains a wide range of ob¬ 
servations for both the Milky Way CMZ and other galactic cen¬ 
tres. In the Milky Way, our model naturally produces the observed 
ring-like stream of gas at ~ 100 pc from the Galactic centre. This 
prediction requires no fine-tuning, and follows simply from the ob¬ 
served shape of the Galactic rotation curve. Second, our model 
naturally explains the broad diversity of star formation rates seen 
in the centres of neraby galaxies, some of which have gas deple¬ 
tion times comparable to those found at larger galactocentric radii, 
and some of which show significantly faster star formation. In our 
model, such diversity occurs naturally as a result of the cycle of ac¬ 
cumulation, starburst, gas clearing, and fading that rapid transport 
produces. It thereby quantifies and greatly expands the scenario of 
episodic star formation proposed by Kruijssen et al. (2014). 

In addition to explaining existing observations, this model 
makes specific predictions for both the morphology and the quan¬ 
titative star formation behaviour of the centres of other galaxies 
when observed at high resolution. While the requisite observations 
have not yet been made, they are well within the reach of ALMA. 
A program to investigate the star formation behaviour of nearby 
galactic centres at high resolution should therefore be a high pri¬ 
ority in the coming ALMA cycles. Such observations would allow 
us to observe galactic centres that are at very different phases of 
the boom and bust cycle than our own relatively quiescent one. A 
sufficiently large sample would even allow us to place statistical 
constraints on the duration of the cycle. This in turn would provide 
strong confirmation that, at least in the most extreme environments, 
star formation is more than a purely local process that proceeds in¬ 
dependent of its galactic environment. 
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